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The law of diffusion governs the transport of heat and charge in most materials and has diverse 
applications ranging from the motion of cells to the prizing of stock options in financial markets. 
Here we investigate the breakdown of diffusion in the transport of fermionic quantum particles on a 
lattice described by a homogeneous Hubbard model. We observe a crossover from diffusive behaviour 
in the center of the cloud to a ballistic motion of atoms in its outer regions. This crossover manifests 
itself in a striking change of the cloud's shape: While it remains round in the diffusive regime, it 
obtains a square shape when the motion becomes ballistic. Surprisingly, the system exhibits a strong 
feedback from the ballistic on the diffusive regions characterized by a universal loss rate of particles 
obeying singular diffusion equations. In addition, the observed dynamics is independent of the sign 
of the interaction, highlighting a novel, dynamic symmetry of the Hubbard model. 



The last years have seen dramatic progress in the 
control of quantum gases in optical lattices [THS] using 
both bosonic and fermionic [4HTQ] gases as well as Bose- 
Fermi mixtures pTHTS] . It has become possible to simu- 
late models of strongly interacting quantum particles for 
which the Hubbard model [14 is probably the most im- 
portant example. A major advantage of these systems 
compared to real solids is the possibility to change all 
relevant parameters in real-time by e.g. varying laser in- 
tensities or magnetic fields. This has led to several stud- 
ies of dynamical properties of both bosonic [T5HT9] and 
fermionic quantum gases JSHTOJ. 

Fermionic atoms have been used to study the equi- 
librium phases of both the attractive [4| and repulsive 
Hubbard model [5Hl] as well as dynamic properties of 
trapped systems [SHTOj. One open, fundamental ques- 
tion in these experiments concerns the timescales needed 
to adiabatically load into the lattice [2QH22] or to achieve 
equilibrium in the lattice [22l|23]. Up to now, all of these 
experiments were performed in the presence of additional 
potentials. 

In this work, it was possible to eliminate all external 
potentials in a two-dimensional (2D) system by compen- 
sating the anticonfining potential of a blue-detuned op- 
tical lattice [6 and to study out-of-equilibrium dynam- 
ics and transport properties in a homogeneous Hubbard 
model. 

To this end, an initially confined atomic cloud was al- 
lowed to expand freely within a homogeneous optical lat- 
tice (Fig. Ilk) for various interactions. Monitoring the 
in-situ density distribution during the expansion allowed 
for several surprising experimental observations: Even 
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FIG. 1. a: Free expansion in a homogeneous lattice. First 
a band- insulator is created in the combination of an optical 
lattice and a strong harmonic trap. Subsequently the har- 
monic confinement is switched off and the cloud expands in 
a homogeneous Hubbard model, b: Observed in-situ density 
distributions. The evolution of the initial density distribution 
(left) crucially depends on the interaction (right). 



small interactions (i) lead to a drastic reduction of the 
expansion velocity of the atomic cloud and (ii) change 
the shape of the expanding cloud. Despite the strong 
effects of the interactions, (iii) we never observe a \/t de- 
pendence of the cloud size, which is expected for simple 
diffusive behaviour, (iv) For strong interactions the core 
width (see below) of the atomic cloud shrinks instead of 
expanding. And, finally, (v) we find that only the mag- 



nitude but not the sign of the interaction matters: the 
dynamics is identical for repulsive and attractive interac- 
tions U. 

In the non-interacting case the atoms expand ballis- 
tically in a continuous quantum walk [24-26 , which is 
also the basis of a recently proposed quantum comput- 
ing algorithm |27j. For long times, this ballistic expan- 
sion results in a square density distribution (see Fig. fTb) 
that reflects the lattice symmetry. Interactions between 
the atoms on the other hand lead to diffusive dynamics 
and allow the particles to regain local thermal equilib- 
rium [28l[29]. As a consequence, sufficiently large atomic 
clouds, where density gradients are small, remain close 
to local thermal equilibrium and their dynamics can be 
described by the laws of hydrodynamics. Thereby the 
cloud retains its initial spherical symmetry (see Fig.flp). 
In the outer regions of these clouds, however, the den- 
sities, and therefore the collision rates, become so small 
that also interacting atoms expand ballistically. 

As diffusion in this system arises from collisions with 
other atoms, the time r between scattering events and 
the diffusion constant D oc r are strongly density de- 
pendent and especially become very large when the den- 
sity of particles n (or holes, 1 — n) becomes very small, 
D{n) (X (n(l — n))~^ (see Supplementary Information 
[30 ). Accordingly, the dynamics is described by a highly 
singular, non-linear diffusion equation: 



dtn{r,t) = V{D{n)Vn) 



(1) 



Such singular (or superfast) diffusion equations have been 
extensively studied in the mathematical literature [31] 
and show many surprising properties. For example, in 
three dimensions (3D), the diffusion equation (Eq.fTj) pre- 
dicts that all particles reach infinity in an infinitesimally 
short time, n{r,t > 0) = if D{n ^ 0) ~ 1/n and if 
lim^^oo n(r, t = 0) =0. This obviously unphysical result 
immediately implies the existence of a strong feedback 
from the ballistic onto the diffusive parts of the atomic 
cloud: The velocity limit imposed by the band structure 
on the ballistic outer regions constrains the maximal ex- 
pansion velocity and thereby holds the diffusive part of 
the cloud together. 

The situation is even more interesting in two dimen- 
sions (2D), where the non-linear diffusion equation pre- 
dicts a universal minimal loss rate of the total number of 
particles N for D{n) = j/n |3Qh32j: 



dtN < -47r7 



(2) 



This loss rate describes the rate of particles reaching in- 
finity and is completely independent of the initial dis- 
tribution. It only depends on boundary conditions at 
r = oo, which can be chosen such that dtN = —Anj [30]- 
[32] . We will show that this loss rate is closely related to 
the rate with which the diffusive part of the cloud emits 
ballistic particles. 

Remarkably, the 2D diffusion equation and the univer- 
sal loss rate have a very interesting geometrical interpre- 
tation ^31^: Consider a curved manifold with a metric 



gij = n{T)6ij. In this case the diffusion equation for 
D{n) = 1/n becomes equivalent to the famous Ricci flow 
'^9ij — ~'^Rij ESI [34:, which describes the evolution 
of a manifold that becomes flatter in time. Here Rij is 
the Ricci tensor describing the curvature of the metric. 
Ricci flows play an important role in mathematics, espe- 
cially for the proof of the celebrated Poincare conjecture 
in three dimensions [34 . Within this interpretation, the 
particle loss rate is universal since dfN is proportional to 
the integral of the curvature over the manifold, which is 
a topological invariant [30] . 



I. EXPERIMENT 



The experiment starts with the preparation of a quan- 
tum degenerate, balanced spin mixture of the two lowest 
hyperfine states \F,mF) = |9/2,-9/2) and |9/2, -7/2) 
of fermionic potassium ^^K in a crossed beam dipole 
trap [6^. Typical clouds consist of N = 2—3 x 10^ 
atoms at an initial temperature of T/Tp = 0.13(2), 
where Tp denotes the Fermi temperature in the harmonic 
trap. Once cooled, the atoms are loaded into a three- 
dimensional blue-detuned optical lattice with lattice con- 
stant A/2 = 369 nm, where their interaction is controlled 
by use of a Feshbach resonance at Bq = 202.1 G [35j. 
The employed loading procedure (see Supplementary In- 
formation [30^) results in a cloud of localized atoms with 
a well-known density distribution, which is independent 
of the interaction between the particles. 

Subsequently, the expansion is initiated by switching 
off the harmonic confinement (see Fig. la). To this end, 
the strength of the dipole trap is reduced by more than 
90%, such that for the horizontal directions the remain- 
ing dipole potential precisely compensates the anticon- 
finement produced by the lattice beams [30 . While the 
vertical motion is expected to be strongly suppressed by 
gravity-induced Bloch oscillations (Oscillation amplitude 
2J/mg < 2 A/2), the atoms are exposed to a homoge- 
neous Hubbard model without additional potentials in 
the horizontal directions. The evolution of the density 
distribution during the following expansion in this quasi 
2D situation was monitored by in-situ imaging along the 
vertical axis of the cloud, thereby integrating over any 
vertical dynamics. Absorption images of the resulting dy- 
namics are shown in Fig. [2] for the case of non-interacting 
particles and in Fig. [4] for various interactions. 

In a second set of experiments, vertical tunneling of 
the atoms during the expansion was suppressed by in- 
creasing the depth of the vertical lattice to 20 Er, with 
recoil energy Er = /i^/(2mA^), thereby realizing several 
layers of independent two-dimensional Hubbard models 
without any influence of gravity. 
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FIG. 2. In-situ absorption images (column density a.u.) of an 
expanding non-interacting cloud in a quasi 2D lattice with lat- 
tice depth 8 Er . The expansion changes the symmetry of the 
cloud from the rotational symmetry of the harmonic trap to 
the square symmetry of the lattice Brillouin zone. At long ex- 
pansion times residual potentials and lattice inhomogeneities 
deform the cloud (see Supplementary Information ^3Qj). 



II. NON-INTERACTING CASE 

In the absence of collisions and additional potentials 

in [30 ) consists 



the Hubbard Hamiltonian (see Eq. Al 



only of the hopping term, which describes the tunneling 
of a particle from one lattice site to a neighbouring site 
with a rate J/h. Each initially localized particle expands 
independently, and the density distribution after an evo- 
lution time t is given by the convolution of the initial 
density distribution with the delocalized probability dis- 
tribution of the individual atoms. During the expansion 
the symmetry of the cloud changes from the rotational 
symmetry of the initial density distribution to a square 
symmetry that is governed by the symmetry of the lattice 
(Fig. |2]). For very long expansion times, this symmetry 
becomes distorted by residual inhomogeneities in the re- 
maining potential and in the lattice depth [30 . 

In order to extract the mean expansion velocity Ve^p 
of the whole cloud we used Gaussian fits [30] to phase- 
contrast images and fitted the evolution of the resulting 
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FIG. 3. Measured cloud size (orange) and deconvolved single- 
particle width (green) of an expanding non-interacting cloud 
in an 8 Er deep lattice. Solid lines denote the quantum- 
mechanical prediction and dashed lines the corresponding 
classical random walk. The inset shows the linear scaling 
of the mean expansion velocity with tunneling J. 



cloud sizes (Fig. [3J orange) by 



m = \Ri 



-^eV 



(3) 
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The deconvolved single-particle width 
^/R{t)'^ — Rq (green) grows linearly in time, as ex- 
pected for a ballistic expansion. In the classical case, 
thermal hopping of a particle (e.g. of a thermalized atom 
on the surface of a crystal) would result in a random 
walk, where at every timestep the particle randomly 
hops to one of the neighbouring sites. The width of the 
resulting density distribution scales as the square root of 
the expansion time. In the quantum case on the other 
hand, a particle tunnels in all directions simultaneously, 
which leads to a fundamentally different dynamics: The 
atoms expand ballistically with a mean squared velocity 
'^exp — (^q)^ which is given by the mean squared group 
velocity Vq = | ^ of all Bloch waves. Here hq denotes 

the quasi-momentum and E{q) = — 2JXli cos(g^|) the 
dispersion relation in the lowest band. For particles that 
are initially localized to single lattice sites, the mean 
expansion velocity is '^exp = V^^f {d: dimension). As 
shown in Fig. [3J the experimentally observed ballistic 
expansion agrees well with the quantum-mechanical 
prediction (solid lines) while a classical random walk 
of the same hopping rate would predict a much slower 
square-root expansion of the single-particle width 
(dashed lines). 



III. INTERACTING CASE 

The ballistic expansion observed for non-interacting 
atoms is in stark contrast to the interacting case, where a 
qualitatively different dynamics is observed: Fig.[4]shows 
in-situ absorption images taken after 25 ms of quasi 2D 
expansion in an 8 Er deep lattice. 
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FIG. 4. In-situ absorption images for various interactions af- 
ter 25 ms expansion in a homogeneous quasi 2D lattice. The 
images show a symmetric crossover from a ballistic expansion 
for non-interacting clouds to an interaction dominated expan- 
sion for both attractive and repulsive interactions. Images are 
averaged over at least five shots and all scales are identical to 
Fig. |2] The bottom line shows the results of a 2D simulation 
of the Boltzmann equation. 



For increasing interaction strengths the center of the 
cloud expands much slower and preserves the initial ro- 
tational symmetry. The observed dynamics gradually 
changes from a purely ballistic expansion in the non- 
interacting case, which results in a square density dis- 
tribution, into a more complex expansion for interacting 
atoms: In the center of the cloud frequent collisions guar- 
antee the system to be close to local equilibrium [28-30^ 
and the diffusion equation (Eq. [l]) applies. As the diffu- 
sion equation is rotationally invariant, the initial spheri- 
cal shape is preserved for most parts of the cloud already 
for moderately strong interactions, see Fig. [T] and |4] For 
a small fraction of atoms in the outer parts of the atomic 
cloud, however, the density is so small that scattering 
becomes rare, so that these atoms expand ballistically. 
Therefore the tails of the cloud show the square symme- 



try characteristic for freely expanding particles (Fig. [4| . 
This initial fraction of ballistically expanding atoms de- 
creases for increasing interaction strengths. During the 
expansion the density gets reduced and, in the limit of in- 
finite expansion times, all atoms are expected to become 
ballistic. We observe the same behaviour irrespective of 
the sign of the interactions. 

The singular diffusion equation (Eq. fl]) would predict 
diverging velocities for atoms at vanishing density. The 
ballistic part of the cloud, however, where the current 
j ^ Vexpn is smaller than the current —D{n)Vn pre- 
dicted by the diffusion equation, obviously limits the 
spread of the diffusive parts of the cloud. This gives rise 
to a strong feedback from the ballistic tails on the diffuse 
core, as discussed in more detail in the Supplementary 
Information [30 . 

For a more quantitative analysis of our results, 
the same experiment was performed in a fully two- 
dimensional (2D) situation, where also numerical simula- 
tions are easier. In this experiment, the vertical lattice is 
kept at a depth of 20 Er during the expansion in an oth- 
erwise identical sequence. In both cases (quasi 2D and 
2D) one observes a spherical, diffusive core surrounded by 
ballistic tails without any qualitative differences between 
the cases. 




FIG. 5. Measured core expansion velocities versus interaction 
for various lattice depths in a 2D situation. The red line 
denotes the result of a numerical calculation (see text). The 
inset shows the doublon dissolution time (see Supplementary 
Information 30 ) and both black lines are guides to the eye. 

The core width Rdf)^ which measures the size of the 
high density core, is extracted from phase-contrast im- 
ages by determining the half width at half maximum of 
the density distribution [30 . Surprisingly, the same sim- 
ple fit function (Eq. [3]) can still describe the time evolu- 
tion (see Fig. 15 in |30 ). The resulting core expansion 
velocities Vc^ which are shown in Fig. |5) decrease dra- 
matically already for interactions much smaller than the 
bandwidth 8 J. This highlights the strong impact of mod- 



erate interactions on mass transport in these systems. 

The core expansion velocities Vc even become negative 
for interactions larger than \U / J\ ^ 3: In this regime, 
the diffusive core dissolves by emitting ballistic particles 
and therefore shrinks in size. Solutions of the diffusion 
equation (Eq. [l]) always show this behaviour for short 
and intermediate times, as the particle current in the 
flanks of the core is always higher than in the center 
of the core, where the density gradient vanishes. The 
slight asymmetry observed at large interactions can be 
attributed to interaction dependent losses due to light- 
assisted collisions during the preparation sequence. 

This pronounced dependence on small interactions en- 
abled us to measure the zero crossing of the scattering 
length {B{a = 0) = 209.1 ± 0.2 G) and to extract a 
new calibration of the width of the Feshbach resonance: 
^ = 7.0±0.2G, see [30. 

For a theoretical description of the expanding clouds 
one needs an approach that can correctly describe both 
the diffusive and the ballistic regime, the probably sim- 
plest one being the Boltzmann equation: 



5t/q + VqVr/q + F(r)Vq/q 



1 



T(n) 



(/q-/q(n)) (4) 



It describes the evolution of the quasi-classical momen- 
tum distribution /q(r, t) as a function of position and 
time in the presence of a force F (see [30^ for details). 
Here the transport scattering time l/r(n), which de- 
scribes the relaxation towards an equilibrium distribu- 
tion /q for given energy (e) and particle densities (n) 
n(r, t) = (n,e), is determined from a microscopic calcu- 
lation of the diffusion constant for small interactions [30] . 
All qualitative features seen in the experiment are well 
reproduced by the numerical results using the Boltzmann 
equation. Fig.[5]shows that the numerical simulations de- 
scribe the drastic collapse of the expansion velocities and 
the shrinking of the core width for strong interactions 
also semi-quantitatively. Quantitative discrepancies be- 
tween experiment and numerics probably arise both be- 
cause the leading order perturbation theory is not valid 
for [/ > J and because the relaxation time approxima- 
tion breaks down in the crossover region from diffusive 
to ballistic behaviour, where the colliding atoms are far 
from thermal equilibrium [30 . 

Even though the core expansion velocity can be quali- 
tatively predicted by a diffusive ansatz, the full quantum 
dynamics is certainly more complex and includes the for- 
mation of entanglement between distant atoms [36]. In 
the case of a sufficiently high initial doublon density the 
free expansion itself could possibly be used to locally cool 
the atoms via quantum distillation processes [37 . 

Our numerical results can be used to study the break- 
down of the diffusion law (Eq. fl]) and the role of the uni- 
versal loss rate (Eq. ^. Fig.p^shows that, for sufficiently 
strong interactions, the diffusion equation can describe 
the behaviour of the cloud center but breaks down in 
the tails. For a quantitative analysis of this effect we 
calculate the cloud size Y^(r^), which, in contrast to the 
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FIG. 6. Calculated density profile multiplied by r^ from a sim- 
ulation of the Boltzmann equation for t = 50 h/ J ^ [/ = 4.15 J . 
Dashed (dash-dotted) lines denote the density profiles along 
the diagonal (the x-direction). The diffusive center is rota- 
tionally invariant and is well described by the diffusion equa- 
tion (blue line). In the ballistic tails, which have a large con- 
tribution to (r^) oc f n{r)r^ dr, the density profiles in diago- 
nal and x-directions differ. The insets show density profiles of 
an expanding atomic cloud (consisting over several indepen- 
dent two-dimensional layers) for various expansion times and 
U = 1.178 J. 



core width, is mostly sensitive to the ballistic parts of the 
cloud (see Fig. [6|. Unfortunately, it cannot be reliably 
extracted from the experimental data due to a combina- 
tion of a small signal-to-noise ratio in the extremely di- 
lute regime and imaging aberrations. The universal loss 
rate suggests that ballistic particles, which are mainly 
responsible for the growth of the cloud size, are created 
at a constant rate. This implies that \/Jr^ grows as t^/^ 
|30j, as confirmed by the numerical simulations (see Fig. 
[1Q| in [30]). Furthermore, even the proportionality con- 
stant can be obtained analytically to within 20 % from 
the universal loss rate of Eq. |2] These numerical and an- 
alytical results strongly support the above picture of an 
emission rate of ballistic particles that is determined by 
the universal loss rate. 

Both for the qualitative and quantitative analysis of 
our results, we assume the relaxation of the system to 
local equilibrium. For very strong attractive or repul- 
sive interactions \U\ ^ J, however, doubly occupied 
sites (doublons) only decay very slowly [23l [381 ES] as 
the missing or excess energy of order U cannot easily be 
transferred to other particles. An important question is 
whether the rate with which the diffusive core melts is 
determined by the decay time of individual doublons or 
whether the latter is fast compared to the former. The 
experimentally observed doublon dissolution times (Fig. 
5| are about an order of magnitude larger than the decay 
time of excess doublons measured recently in a half filled 
situation [23 in 3D at comparable interactions. Fur- 
thermore, they match very well the typical timescales 
for melting of the diffusive core of the cloud observed 
within our simulations. We therefore conclude that for 



i 



the parameters used in our experiment the doubly oc- 
cupied sites remain in local equilibrium in the diffusive 
regime. 

Surprisingly, we measure identical density profiles and 
expansion rates for repulsive and attractive interactions 
of the same strength (see Fig. [s] and |4|, although one 
would expect that repulsive (attractive) interactions lead 
to a positive (negative) pressure and therefore an in- 
creased (reduced) expansion rate. While this simple pic- 
ture indeed applies for atoms with a quadratic dispersion 
p^ /2m^ it is not correct for the Hubbard model with a 
periodic dispersion relation £^(q). In combination with 
an initial state that consists of a homogeneously filled 
Brillouin zone, where quasi-momenta with positive effec- 
tive mass {d'^E{q)/dq^ > 0) and negative mass (hole- 
like states) are equally occupied, the high symmetry of 
the kinetic energy guarantees identical density profiles for 
interactions U and —U: As both initial state (localized 
atoms) and observable (density) are invariant under time 
reversal symmetry, the density dynamics is unchanged 
by the transformation H ^ —H. In addition, they also 
remain invariant under a boost q -^ q + (7r,7r,7r) x -|, 
which changes the sign of the hopping J ^ —J. In com- 
bination, these two invariances give rise to the observed 
dynamical symmetry (see also [30]). 



IV. CONCLUSION 

Ultracold fermions in optical lattices offer many novel 
possibilities to study non-equilibrium dynamics as they 
allow for a full real-time control of most relevant param- 
eters. Particularly they can implement so-called quan- 
tum quenches, where the Hamiltonian of the system is 
changed instantaneously. We studied the expansion of a 
cloud of initially localized atoms in a homogeneous Hub- 
bard model and observed the crossover from a ballistic 
expansion at small densities or vanishing interactions to 



hydrodynamic expansion in the interacting case. The 
feedback between the diffusive and ballistic parts of the 
cloud controls the expansion: the diffusive core slowly 
emits ballistic particles which in turn hold the diffusive 
part of the cloud together and regularize the otherwise 
singular diffusion in the tails. We observed identical be- 
haviour for both attractive and repulsive interactions, 
highlighting the high symmetry of the kinetic energy in 
the Hubbard model. 

The surprisingly large timescales of mass transport in 
an interacting Hubbard model set lower limits on the 
timescales needed both to adiabatically load the atoms 
into the lattice and to cool the system in the lattice [40] . 
They are therefore of paramount importance for all at- 
tempts to create complex, strongly correlated many-body 
states like Neel-ordered states in these systems. The 
method of directly measuring the expansion velocity can 
be generalized in a straightforward way to more complex 
quantum states including metallic and Mott-insulating 
states in the repulsive Hubbard model [6 or the pseu- 
dogap regime in the attractive Hubbard model [4 . In 
addition, the effects of various disorder potentials on the 
two-dimensional dynamics can be studied. 
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Appendix A: Supplementary Information 

1. Experimental sequence 

The experiment starts with the preparation of a band- 
insulating state of a balanced spin mixture of fermionic 
potassium ^^K, as described in previous work [6]: 
Through evaporative cooling in a red-detuned crossed 
dipole trap a quantum degenerate mixture of the two 
lowest hyperfine states of potassium was reached with 



atom numbers of N = 1 — 1.5 x 10^ atoms per spin state 
at a temperature of T/Tp = 0.13(2), where Tp denotes 
the Fermi temperature in the harmonic trap. The trap- 
ping frequencies of the dipole trap were then increased 
to 27r X 100 Hz (27r x 400 Hz) in the horizontal (vertical) 
directions. 

Subsequently, a simple-cubic blue-detuned 3D optical 
lattice (A = 738 nm) is ramped up linearly to a depth 
of 8Er {lEr = K^ /{2m\^)) in 56ms while the magnetic 
field is held at 209.1 G, which corresponds to vanishing 
interactions (see A 6). This loading procedure results in a 



band-insulating state surrounded by a thin metallic shell 
at a compression oi Et/12J = 1.8, using units and con- 
ventions of Ref. [6 . In the next step, the tunneling rate 
J is reduced to J = /i x 23 Hz by linearly increasing 
the lattice depth to 20 E^ in 200 /is, a timescale that is 
slow enough to avoid excitations into excited bands, but 
fast compared to tunneling within the lowest band. Due 
to this reduced tunneling rate the density distribution 
is essentially frozen out during the following 40 ms mag- 
netic field ramp {Bdyn = 206 — 260 G), which sets the 
interaction for the expansion. Combined with the strong 
harmonic confinement, this leads to a dephasing between 
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FIG. 7. Experimental Sequence. Starting with a degenerate 
Fermi gas in the dipole trap a non-interacting band insula- 
tor is created. During a freeze-out period the atoms localize 
to individual lattice sites and the desired interaction is set 
without altering the density distribution. Subsequently the 
harmonic confinement is switched off and the cloud expands 
in a homogeneous Hubbard model. 



different lattice sites and effectively localizes all particles 
to individual sites [41 . 

In total, this sequence produces a cloud of localized atoms 
with a well-known density distribution that is indepen- 
dent of the interaction between the particles. The expan- 
sion is initiated by lowering the lattice depth in 200/is 
to values between AEr and 15 Er while simultaneously 
switching off the harmonic confinement. To this means, 
the strength of the dipole trap is reduced by more than 
90%, such that the remaining potential precisely compen- 
sates the anticonfinement produced by the lattice beams 
(see Sec.[A4). 



2. Theoretical Analysis: Equilibration, Superfast 
Diffusion and the BreakdoAvn of Hydrodynamics 

The hydrodynamic approach is a powerful method to 
describe the dynamics of interacting liquids. It can be 
used in many different situation to model the properties 
of gases (e.g. for the weather forecast), liquids or plasmas 
both in the classical and quantum regime. As we will 
show below, in our system the hydrodynamic approach 
is essential to understand the properties of the center 
of the interacting atomic cloud where the dynamics is 
diffusive. It breaks down, however, in its tails where 
scattering processes of the diluted quantum gas are rare 
and the dynamics is ballistic. Surprisingly, there is a 
strong feedback from the ballistic tails to the diffusive 
core. This feedback can be traced back to a peculiar 
type of superfast diffusion. 

In the following we will first show that a simple hy- 
drodynamic approach based on a diffusion equation is 
expected to appropriately describe the dynamics in the 
center, but not close to the border of the atomic cloud, 
where densities are very low. Due to the specific form of 
the nonlinear diffusion equation, we will argue that even 
in the center of the cloud the properties are strongly in- 



fluenced by the ballistic boundary. We will therefore ana- 
lyze and describe the interplay between the diffusive core 
and the ballistic tails using a Boltzmann equation. 



a. Diffusive Regime 

The hydrodynamic approach is mainly based on the 
fact that scattering processes tend to equilibrate the sys- 
tem locally. This equilibration occurs on a time scale 
'?"scatt7 the time between collisions of the particles, which 
depends on the microscopic details of the scattering pro- 
cesses. Similarly, there is an associated length scale 
^scatt 5 which can be identified with the typical distance an 
atom travels between two collisions. If the system now 
is probed on time and length scales much larger than 
Tscatt and ^scatt 7 the microscopic details of the scattering 
processes are not important any more and the system 
is locally close to equilibrium. The local equilibrium de- 
pends only on a small number of conserved quantities 
that remain invariant under scattering. 

Starting point of our analysis is the Hubbard Hamil- 
tonian 



H 



-^E 



^ia^jcr 



u 






rii^riii 



(Al) 



(iJ) 



where q^ creates a fermion with spin a =t/i on lattice 
site i and where riia = cJo-^^cr- The first term describes the 
tunneling with rate J from one site of the square lattice 
to a neighbouring site, while the on-site interaction of the 
atoms is described by U. 

In this section we use units where h = 1 and where the 
lattice constant A/2 is set to 1. This has the advantage 
that the density per area and the density per lattice site 
take the same values. 

For atoms in an optical lattice, momentum can be 
transferred to the optical lattice by a combination of 
particle-particle and Bragg scattering and there remain 
only three conserved quantities in the Hubbard model: 
the particle number of each of the two species and the 
energy. As the local magnetization always vanishes in 
our setup, the hydrodynamic equations can be formu- 
lated in terms of two densities only: The number density 
n(r,t) and the energy density e(r, t). In the following, 
we combine both quantities into a vector n = (n, e) for 
notational convenience. The hydrodynamic equations are 
simply derived from the continuity equation Q^n+Vj = 0. 
The energy and particle currents j vanish in thermal equi- 
librium where n = const. For smooth variations of n one 
can expand the current j to linear order in Vn to obtain 
the diffusion equation 



dtn = V(D(n)Vn) 



(A2) 



where D is a two-by-two matrix, describing the energy- 
and particle diffusion constants and cross terms which 
arise as energy gradients induce particle flow and vice 
versa. D(n) can only be obtained from a microscopic 



calculation (see below). For the following discussion it 
will be essential that for low densities D{n) ex 1/n, as the 
scattering rate is proportional to the number of scattering 
partners (see Sec. A 2c). 

Note that, in one dimension, the hydrodynamic ap- 



b. Breakdown of Hydrodynamics 



proach cannot be used for the Hubbard model (Eq. Al), 



since there exists an infinite number of further conserva- 
tion laws that strongly restrict the possibilities for equi- 
libration. 



The diffusion equation (Eq. A2 ) is valid under the con- 
dition that the scattering rates are large compared to 
the expansion rate 1/rexp.. A quantitative version of this 
condition can be obtained by linearizing the Boltzmann 
equation, see Sec. |A2c[ around the local equilibrium so- 
lution /q ^ /^(n) + (5/q with J/q (X Tgcatt • To leading 
order in Tgcatt one obtains the diffusion equation for n. 
By comparing the subleading correction to the current to 
the leading terms, one obtains the condition 



1 1 
> 

'^scatt '^exj 



|V(D(V(DVn)))| 
IDVnl 



(A3) 



for the validity of the diffusion equation (Eq. A2 ) . Ap- 
proximating gradients by the inverse radius of the cloud, 
1/R ^ 1/A^^/^, and using D ^ '^^Tgcatt where '^ '^ J is 
a typical velocity, N the total number of particles in the 
experiment, and J the hopping rate, one obtains that the 
diffusion equation is valid for 



1 
> 

'^scatt 

n> 



J 



iVl/3 
[/2JV1/3 



(A4) 
(A5) 



The second equation is valid for small interactions U 
only and makes use of t he low density approximation 

l/'^scatt ^ nU'^/J^ see Sec. A 2c 



An alternative, less rigorous approach to estimate the 
validity range of the diffusion equation assumes it to be 
valid as long as the currents jbaii ^ vn for a ballisti- 
cally expanding cloud are much larger than those ob- 
tained from the diffusion equation. 



jdiff = |I^(n)Vn| < jbaii 



This results in 



^scatt 



> w- 



iVnl 



(A6) 



(A7) 



consistent with Eq. A3 and Eq. |A4[ 

In the thermodynamic limit (defined as usual by A^ ^ 
oo while scaling the initial strength of parabolic confining 
potentials with 1/A^^/^) the diffusion equation appears to 
be valid for the whole cloud according to Eq. |A5| In the 
experimental case of 10^ particles and moderate interac- 
tions it should be possible to describe almost all atoms 
by the diffusion equation Eq. |A2| Surprisingly, this con- 



clusion turns out to be incorrect, as can already be seen 
by analyzing the solutions of the diffusion equation. 



According to Eq. |A3| and the estimate Eq. |A5[ the 
diffusion equation is not valid in the tails of the cloud 
where, due to an extremely small scattering probabil- 
ity, a crossover to ballistic motion is expected. It is, 
nevertheless, very instructive to study the properties of 
the diffusion equation Eq. |A2[ since, as discussed above, 
the largest part of the cloud remains diffusive. Here we 
are mainly interested in the low-density regime, where 
D{n) ex 1/n, since we want to analyze the breakdown of 
diffusion in the tails. Heat currents are negligibly small 
for our experiment. Also the dependence of the diffusion 
constant on energy can be neglected as the local tem- 
peratures are always much larger than the bandwidth, 
\T/J\ ^ 1, as we have checked numerically. Therefore, 
we focus on the simple diffusion equation 



dtn{r,t) 



1 



n{r,ty 



:Vn{r,t) 



(A8) 



with (5 = 1 in the experiment. A prefactor parameterizing 
the strength of interactions can be absorbed in a rescaled 
time variable. 

For 5 > 0, Eq. |A8| is often called fast diffusion equation, 
for (5 > 1 superfast diffusion equation [31^, and for J = 1 
also the name logarithmic diffusion is used as Vlogn = 
^Vn. 



Nonlinear diffusion equations of the form of Eq. A8 



investigated intensively in the mathematical literature, 
for a review see ^Slj. Also in physics ^42ti46 diffusion 
laws with 5 = 1 have been discussed in a small number 
of publications to describe e.g. the effect of magnetic mir- 
roring in plasmas [44] , the kinetics of ion exchange j3Sl , 
or the spreading of liquid drops on solids in a regime dom- 
inated by van-der-Waals forces [46 . A clear picture on 
how singular hydrodynamic equations and the associated 
unphysical results have to be regularized and interpreted 
appears to be missing both in the physical and the math- 
ematical literature. Good examples for the associated 
theoretical and experimental challenges can, for exam- 
ple, be found in the literature on the dynamics of liquid 
drops on surfaces. There, mesoscopic precursor films and 
single-molecular layers apparently play an important role 
for the regularization of hydrodynamic singularities but 
their infiuence on the dynamics remains unclear [47 . 

Here we will describe some of the basic results relevant 
for our discussion, mainly based on a review by Vazquez 
J3T] . A lot of insight into the solutions of fast diffusion 
equations can be obtained by studying simple scaling so- 
lutions (the celebrated Barenblatt solutions [48,) of the 
form 



n(r,t) = ^/(r/t") 



(A9) 

Inserting this ansatz into Eq. |A8| one finds in d dimen- 
sions 

^ ^ {5 = 1) (AlO) 



2-d5 2-d 



10 



indicating that for d = 3 no such scahng solution exists 
while d = 2 constitutes the marginal case. 

It is also instructive to check the conservation of the to- 
tal number of particles by integrating the diffusion equa- 
tion Eg. |A8| from r = to r = rmax -^ oo: 



^d-l 



dfN (X 



-n\r) 



(All) 



One immediately realizes that the surface term does 
not vanish in general, when n ^ at the edge of 
the cloud: particles are lost at r = oo. The impor- 
tance of the boundary term obviously depends both on 
S and d. For our case, (5 = 1, a vanishing surface term 
{dlnn/dr < —c/r^~^ for all c > 0) is only possible for 
d < 2 when J l/r^~^ diverges. A loss of particles is 
therefore unavoidable for d > 2^ S = 1. 

This implies several quite dramatic consequences 
which have been worked out in the mathematical litera- 
ture [3TJ |32l |49] . As a first surprising result, for a given 
initial condition n{r^t = 0) with n{r ^ oo,0) = 0, the 
diffusion equation Eq. |A8| does not even have a unique 
solution for ^ = 1! Due to the finite surface contribution 
one can specify an arbitrary time-dependent current at 
infinity. For d = 1^ S = 1 (or, more precisely, d < 2) it 
is, however, possible to find solutions with a conserved 
number of particles, e.g. n{r,t) = 2t/(r^ + ('^^)^) with 
arbitrary velocity v, consistent with a = 1 obtained from 
Eq. |A1Q[ In d = 3, in contrast, a is negative indicat- 
ing the nonexistence of a scaling solution. Indeed, it 
has been shown [STJ 09] that for arbitrary initial con- 
ditions with n{r -^ oo,0) = a solution of the diffusion 
equation does not exist for 5 = 1: all particles reach 
r = oo in an infinitesimally short time. The situation 
is most interesting in two dimensions {d = 2), where 



a ^ oo according to Eq. AIQ This case can easily be an- 
alyzed for radially symmetric initial conditions by rewrit- 
ing n(r, t) = n(ln(r),t)/r^. For ^ = 1 one finds h{u^t) to 
follow the one-dimensional fast diffusion equation with 
^ = 1. The total number of particles is not conserved 
but instead drops in time, 

dtN = -47r + 27r lim ^^f^ < -47r lim [nD{n)] (A12) 

where we have recovered the correct prefactor in the last 
term. Interestingly, there is a universal minimal loss rate 
which is independent of the size and shape of the cloud 
[3TJ[32]. This is due to h\u -^ oo) < 0, because n{r^t) 
has to decay faster than 1/r^ to obtain a finite particle 
number 27t J n{r^t)rdr. The minimal loss rate can, for 
example, be obtained by solving the diffusion equation for 
an initial condition ne{r^ t = 0) = no{r) -\- e with no(r -^ 
oo) = 0. For finite e the particle number is conserved as 
D{n) is bounded from above. The limit liiRe^o ne{r^t)^ 
however, describes a solution where the minimal loss rate 
47r is realized [31]. After a finite time tmax proportional 
to the total number of particles in the two-dimensional 
sheet, tmax ^ U'^ / JN2d^ all particles have disappeared. 



The diffusion equation predicts them all to have reached 
r = oo. 

In addition, the factor 47r exhibits a direct geometrical 
interpretation using the relation to the Ricci flow dis- 
cussed in the main text: For the two-dimensional metric 
Qij = n8ij the total particle number J n = J ^/\detg\ is 
given by the area of the manifold M, and the particle 
loss rate from a Ricci flow is given by an integral over 
the curvature K, dfN = —2jj^K. For a compact mani- 
fold without boundaries, the integral over the curvature, 
Jj^ K^ is completely fixed by topology. According to the 
Gauss-Bonnet theorem it is given by 27r times the Eu- 
ler characteristic, which is 1 in our case, leading to the 
rate —Att. The extra surface term in Eq. A12 also has 
a similar interpretation in terms of an integral over the 
geodesic curvature over the boundary of the manifold. 

The mapping for the two-dimensional to the one- 
dimensional problem also shows that the exponent a = 
oo obtained from Eq. AIQ for ^ = l,(i = 2 can be in- 
terpreted as an exponentially fast growth of the cloud, 
r ~ exp(7;t). 

Obviously, the results obtained from the analysis of 
the diffusion equation are "unphysical" in the sense that 
in reality there is no particle loss and the cloud cannot 
expand faster than ballistically in the long-time limit. 
Nevertheless, the analysis of the diffusion equation shows 
that there is an unexpected feedback mechanism between 
diffusive and ballistic regime: The sensitivity of the diffu- 
sion equation to boundary conditions clearly shows that 
the ballistic tails also control the dynamics of the dif- 
fusive core. In three dimensions the diffusive core is 
literally held together by the ballistic tails as the dif- 
fusion equation alone would predict that all particles 
vanish in an infinitesimally short time. Similarly, in 
the two-dimensional case, where the diffusion equation 
predicts a universal loss rate of particles at infinity, the 
crossover from diffusive to ballistic behaviour will affect 
both regimes. 

From a more general point of view, the experimentally 
realized case, d = 2^5 = 1 constitutes a multicritical 
point: For 5 = 1^ d = 2 represents the critical dimen- 
sion, because particle loss is unavoidable for all d > 2, 
as shown above. But also in higher dimensions 5 = 1 
defines a special point, as can already be seen by using 
e.g., an exponential density profile n ^ e~^^ in Eq. A12 



which yields a diverging (vanishing) particle loss for 6 > 1 
{S < 1). Indeed, for (5 > 1 and d > 2 all particles vanish 
instantaneously at infinity, while for 1 > S > 2/d solu- 
tions exist for finite time intervals ^Slj. It is therefore 
both experimentally and theoretically an interesting goal 
for the future to investigate the interplay of ballistic and 
diffusive motion for d = 3, 5 = 1. 



c. Boltzmann equation 

In the preceding section we have shown that the feed- 
back from the ballistic tails on the diffusive core of the 
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atomic cloud controls its expansion. Consequently, an 
appropriate theoretical description requires a model that 
can describe both the diffusive and the ballistic regime 
as well as the crossover between the two. These require- 
ments are met by the Boltzmann equation. A full solu- 
tion of the space-and time-dependent Boltzmann equa- 
tion is, however, numerically very demanding. We there- 
fore use a simplified version, which is constructed to be 
exact both in the diffusive and the ballistic regime and 
is able to provide a decent qualitative (but not quanti- 
tative) description of the crossover regime: We employ 
a version of the relaxation time approximation that con- 
serves both energy and momentum. Boltzmann equa- 
tions [50] describe the evolution of the quasiclassical mo- 
mentum distribution /q(r, t) as a function of position and 
time. Here, the equation is given by 



9t/q + VqVr/q + F(r)Vq/, 



1 



qjq 



r(n) 



(/q-/q(n)) 



(A13) 



where Vn 



VqCq 



denotes the velocity of the particles 



and F(r) describes a force term arising from the inter- 
action of the atoms. For the initial condition we take 
into account that coherences between neighboring sites 
are efficiently destroyed during the preparation of the 
initial state, as described in Sec. |A 1[ Therefore we use 
/q(r,t = 0) = n(r) where n(r) is the density distribution 
of non-interacting fermions with an entropy per particle 
of 1.141 kB- This entropy corresponds to an initial tem- 
perature oiT/Tp ~ 0.125 in the harmonic trap, which is 
compatible with the experiment. In the force term, we 
use the Hartree approximation, F = —UVn^ where n is 
the density per spin, but our numerical results show that 
this term is relatively small. In order to obtain the cor- 
rect hydrodynamics, the collision term on the right-hand 
side of Eq. |A13| needs to conserve both particle number 
and energy, and to correctly describe the relaxation to- 
wards equilibrium. This is achieved by setting 

/°(n) = l/(exp[(eq - ^^(r, t))/kBT{v, t)] + 1) (A14) 

where /J^{r^t) and T{T,t) are chosen such that both 
E, /q - /^(n) = and E, ^q (A - /^(n)) = 0. The 
effective temperatures T(r) obtained from this procedure 
are, for our initial conditions, at least an order of magni- 
tude larger than the bandwidth. This implies that energy 
conservation and the diffusion of energy only play a mi- 
nor role in the experiment, as we checked numerically. 

The effective scattering rate l/r(n) should be chosen 
such that the correct diffusion constant D{n) is obtained. 
Unfortunately, we are not aware of any numerical or an- 
alytical method to calculate D{n) = r(n)(vq^)/(i for 
the two-dimensional Hubbard model, even in the limit 
T ^ oo, where all thermodynamic properties are ex- 
actly known. Note that, for example, dynamical mean 
field theory, which successfully describes thermodynamic 
properties [6], should not be used to calculate transport 
properties quantitatively, as it neglects vertex corrections 
that are quantitatively important even for T ^ oo. 




FIG. 8. Transport scattering rates as a function of den- 
sity for different temperatures. Note that the transport 
scattering rate (in contrast to the single-particle scattering 
rate) becomes exponentially small in the low-density, low- 
temperature regime where Umklapp scattering is suppressed. 



We therefore calculate the diffusion constant pertur- 
batively to leading order in the interaction. This can be 
done by computing the conductivity a from the transla- 
tionally invariant Boltzmann equation [50] , using the full 
collision term with golden-rule transition rates. Subse- 
quently, one can obtain D{n) and r(n) from the Einstein 
relation D = ad/jj/dn. 

In order to avoid the solution of a complicated inte- 
gral equation, we employ a variational solution of the 
Boltzmann equation (see [S^ for details) using /q = 

/q ~ ^/q/^^q Xli <^iCq with four Variational parameters 



and 



.(1) 



.(2) 



.(3) 



.(4) 



L-q — '^q? '-q — '^q'^q? '^q — Qx") '-q 

{ir/a — Qx) mod 27r/a. Here c^^"^ and c^^^ are chosen to 
enable the calculation of both thermal and charge diffu- 
sion constant as well as the cross terms. The remaining 
terms, c^^^ and c^^\ on the other hand, are essential to 
correctly describe the low temperature limit: For a low 
density of particles or holes, the conductivity, and there- 
fore the diffusion constant, grows exponentially for low 
T due to an exponential suppression of Umklapp scat- 
tering processes [50 . Since we determine l/r(n) such 
that we recover the correct diffusion constant, this effect 
is fully included. Note that it is not captured by a more 
conventional version of the relaxation time approxima- 
tion [50 , which neglects vertex correction and identifies 
l/r(n) with a single-particle relaxation rate instead of a 
transport scattering rate. Although the variational ap- 
proach only gives a lower bound for the diffusion con- 
stant, we expect it to be accurate within a few percent 
for small interactions, as we have checked by reducing 
the number of variational parameters. The density de- 
pendence of the resulting scattering rate is shown in Fig. 
|8] for various temperatures. For £^ = or, equivalent ly, 
T ^ oo we obtain l/T{n,E = 0) ^ 0.609 n(l - n)lPj J . 
To solve the resulting Boltzmann equation (Eq. A13), a 
simple Runge-Kutta integration is used, discretizing both 
the position and momentum variable. 
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FIG. 9. Cloud size Y^(r^), obtained from a simulation of the 
Boltzmann equation, as a function of expansion time for var- 
ious interaction strengths. The initial conditions are chosen 
to match the experimental parameters. 



For all quantitative comparisons with the experimental 
results we also take into account the necessity to aver- 
age over a stack of independent two-dimensional systems 
with different atom numbers. This averaging, however, 
changes the result only slightly since the main contribu- 
tions arise from the central layers. 



d. Interplay of Diffusive and Ballistic Regime 

The analysis of the diffusion equation reveals that, in 
order to appropriately describe the diffusive core of the 
atomic cloud, one has to include the effect of the ballistic 
outer regions as well. As the Boltzmann equation de- 
scribes both regimes and their crossover, it can be used 
to determine the characteristics of their mutual influence. 

Both in the experiment and in the simulation the diffu- 
sive and ballistic parts of the atomic cloud can easily be 
distinguished. The diffusion equation is rotationally in- 
variant, while, in contrast, the cloud becomes quadratic 
in the ballistic regime due to the square symmetry of 
both dispersion and Brillouin zone, see Fig. 1 and 6 of 
the main text. A direct comparison of the density pro- 
files obtained in the simulations and in our experiment 
demonstrates a very good agreement on all qualitative 
features. Using scattering rates that were obtained to 
leading order in perturbation theory, a full quantitative 
agreement could not be reached. This is discussed in 
detail in Sec. lA 2 el below. 

Diffusive and ballistic regimes couple to each other in 
at least two ways. On the one hand, the diffusive core of 
the cloud "feeds" the ballistic tails: the diffusive dynam- 
ics determines how many particles emerge in the low- 
density tails of the atomic cloud. On the other hand, 
the diffusive center is held together by the ballistic tails 
of the cloud, since the cloud would expand much more 
rapidly within the diffusion approximation alone. One 
way of perceiving the interplay of the two regimes is to 



study different definitions of the radius of the cloud. Fig. 
|6] of the main text shows that the mean-square displace- 
ment (r^) is mainly determined by the ballistic tails of 
the cloud. In contrast, for not too weak interactions, the 
core width (half width at half maximum HWHM) lies 
in the diffusive core of the cloud. As the diffusive core 
melts and loses particles to the ballistic outer regions it 
is therefore not surprising that, for large interactions, the 
core width shrinks as a function of time (both in theory 
and experiments), (r^), on the other hand, always grows 
(see Fig.|9|. 

An interesting question concerns the long-time asymp- 
tote of the cloud radius: Does it show power-law be- 
haviour r ~ t" in at least some regime? A related ques- 
tion is whether the universal loss rate of particles ob- 
tained from the hydrodynamic theory in d = 2 has any 
direct consequences. 

To study these questions, it is useful to first investigate 
the behaviour of the ballistic part of the cloud. The dif- 
fusive core can be viewed as a source of ballistic particles. 
In order to obtain a qualitative impression of the shape 
and effective radius of the ballistic part of the atomic 
cloud, we can use a simple model: At r = 0, ballistic 
particles are created with rate W{t) and momentum dis- 
tribution nq, which is assumed to be time independent. 
A ballistic particle created at time to at r = with ve- 
locity Vq can be found at later times t at r = Vq(t — to). 
Thus, we obtain the following density distribution of the 
ballistic atoms: 



f d^q f 



W{to)S^{r- 



■ vq(t - to)) 
(A15) 



More realistically, one can assume that the particles are 
created at a finite radius Tc, which avoids, e.g., an un- 
physical logarithmic singularity on the diagonal x = y. 
Since we are mainly interested in the limit r ^ Tc, we 
ignore this effect. As shown in Fig. [6J the mean radius 
Y^(r^) is dominated by the ballistic part of the cloud. 
Therefore, for a single two-dimensional sheet with N2d 
atoms per spin, we can estimate 



{r% - {r%=o 



/r2nbaii(r,t)rf2r 



(v^) 



N2d 

!oW{to){t-tofdto 
J^ W{to)dto 



(A16) 



Finally, for i — > oo, all particles become ballistic, thus 

/o°° W{to)dto = N2d. 

Since the diffusion equation predicts the disappearance 
of particles at infinity with a universal, time-independent 
rate, we define a time t^ax such that W{t) ^ A^2d/^max 
for t < tmax and W{t) = for t > tmax- This leads to: 



(r^), - {r2),=o « K) 



2v I iV(3Wx) iort<t 



for t > tmax 

(A17) 
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FIG. 10. Cloud size y^(r^) as a function of t^^^ (in rescaled 
units). While the core width (half width at half maximum 
HWHM) of the atomic cloud measures the radius of the dif- 
fusive core, the average yjr^ is mainly sensitive to the bal- 
listic tails. From the numerical solution of the Boltzmann 
equation, we find that \/ {r'^) grows as t^''^ for long times, 
see [30]. Remarkably, the prefactor of the t^''^ law (see inset) 
can be calculated analytically (dashed line) from the univer- 
sal loss rate of the diffusion equation discussed in Eq. |A12| 
For weak interactions and much longer times, t ^ tmax, the 
whole cloud is ballistic with yjjr^ oc t. 



Interestingly, precisely this behaviour is observed in 
our numerical solution of the Boltzmann equation. In 



Fig. 

Vtm 



10 



the size of the cloud, r = y (r^), scaled by 
is plotted as a function of t^/^. Here, as be- 
fore, V = Y^(vq^) denotes a typical ballistic velocity. Re- 
markably, one does not only obtain r/{vtmax) ~ t^^'^ for 
t < ^max but also that the prefactor is independent of the 
interaction strength! 



Eq. A12 yields an estimate of the prefactor of the t^/^ 
behaviour 



N^ 



2d/{^7r\im[nD{n)]). 



(A18) 



This prefactor differs by less than 20 % from the one ob- 
served in the numerical solution of the Boltzmann equa- 
tion in the high-density limit, see Fig. IlO) A similar be- 



haviour is also found for initial conditions with a much 
smaller central density, n <C 1, as long as t < tmax (not 
shown) . 

Overall, the agreement of the numerical simulations 
and the simple model (Eq. A17) is rather surprising as 



the diffusive behaviour is strongly modified by the ballis- 
tic outer regions, which prohibit a too fast expansion of 
the cloud. Here one should note that from the present nu- 
merical data we cannot rule out an additional very weak 
logarithmic dependence of the prefactor on cloud size, 
interaction strength, or time, which would not be sur- 
prising as d = 2 is the 'critical' dimension of the problem 
as discussed above. 



e. Origin of quantitative discrepancies of theory and 
experiment 

While the numerical results qualitatively reproduce all 
experimentally observed features, a quantitative compar- 
ison shows that the experimentally observed expansion 
rates of the cloud are consistently larger than the esti- 
mates obtained from our numerical results. Note that 
there are no fitting parameters involved as both the hop- 
ping rates and the interaction strengths are known. The 
error bars in J and U cannot explain the observed dis- 
crepancies. 

There are several candidates that could explain the dis- 
crepancies between theory and experiment: In the non- 
interacting case, where there are no theoretical uncertain- 
ties, the agreement of theory and experiment concerning 
the size of the cloud is very good and pronounced discrep- 
ancies in the shape only occur for large clouds, where 
residual potentials and inhomogeneities in the hopping 
rates become important, see Fig. [2] in the main text and 
Sec. |A4| below. Due to the much smaller core widths, 



however, we do not believe that these effects are impor- 
tant in the interacting case. 

A second problem might be the variational estimate 
for the diffusion constant, which gives only a lower bound 
for D. While the experimentally observed cloud sizes are 
indeed larger than the theoretical predictions, we do not 
believe that errors due to the bound are able t o explain 
the discrepancy (see discussion in section A2c above). 



A more serious problem is the use of perturbative for- 
mulas for the diffusion constant, which are only valid 
for weak interactions, formally U <C 4J. A quantitative 
discrepancy for large interactions is therefore not surpris- 
ing. We were, however, not able to identify a clear trend 
whether the theory fits better for weaker interactions. 
Consequently, this argument does not give a convincing 
explanation either. 

In our opinion, the most likely origin of the observed 
discrepancy lies in a breakdown of the relaxation time ap- 
proximation in the crossover regime between the diffusive 
core of the cloud and its ballistic tails. Deep in the bal- 
listic regime, the scattering rate is negligibly small. Sim- 
ilarly, deep in the diffusive regime, the relaxation time 
approximation should be valid: To good approximation, 
a local equilibrium is established and we have determined 
l/r(n) in such a way that the correct diffusion constants 
are reproduced. Yet in the crossover regime, the distri- 
bution functions /q are far from equilibrium. Here, the 
relaxation time approximation should result in an error 
of order one. While this crossover regime only affects a 
small portion of the cloud, we know from our analysis of 
the diffusion equation in Sec. |A2b| that boundary condi- 
tions are very important for diffusion equations. Conse- 
quently, the crossover regime will have a strong influence 
even on the core width of the cloud, which, for moderate 
interactions, is determined from the portion of the cloud 
where the diffusion equation should be valid. 
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3. Dynamical U 



-U symmetry of the Hubbard 
model 



In the fohowing we address the observed dynamical 
symmetry between repulsive and attractive interactions 
in the fermionic Hubbard model. As discussed in the 
main text, this symmetry arises from the high symmetry 
of the tight binding dispersion ^q = — 2JX]i cos(g'i|). 
In order to understand this, let us consider a single mo- 
mentum component, say hqx^ whose dispersion and group 
velocity distribution are plotted in Fig. pT| 

The maximum group velocity 'Umax = ^^^(|; a"^) ^^^" 
responds to the state with quasi-momentum hq^ — fm/X 
and energy £^g = 0, while states with finite energy, either 
positive or negative, have a smaller velocity. Using this 
group velocity distribution, the symmetry of expansion 
under the change of the sign of U can intuitively be un- 
derstood as follows: Attractive interactions reduce the ki- 
netic energy of the particles, which slows them down. Re- 
pulsive interactions increase their kinetic energy which, 
surprisingly, slows them down as well. 




-Vmax 



FIG. 11. Single particle energy Eq and corresponding velocity 
Vq in the Hubbard model. 



In the following we will turn this qualitative argument 
into a rigorous theorem revealing this counter- intuitive 
dynamical symmetry of the Hubbard model: 

We consider a coherent dynamical evolution arising 
from two Hubbard-type Hamiltonians that differ only in 
the sign of the interaction term: 



H± 






^icT^Jcr 



±UY,^^^^^i (A19) 



In addition, we define two symmetry operations: The 
time reversal operator Rt and the 7r-boost operator Bq^ 
which translates all quasi-momenta by Q = (tt, tt, tt) x I. 
Applying the time reversal operator turns the wavefunc- 
tion into its complex conjugate [51 or, equivalently, mod- 
ifies the time evolution operator: 

Rte-'^^R\ = e'^^ (A20) 

The action of the 7r-boost operator is given in second 



quantized notations by 

BqCpBQ = Cp^Q (A21) 

or, in real space: 

BQCjBQ = e'^'^^Cj (A22) 

We now formulate the following general theorem: 

// both the initial state |^o) ^^^ the experimentally 
measured quantity O are invariant under both time re- 
versal and 7T-boost, the observed time evolutions 



are identical: (0(t))+ = (0(t))_. 



(A23) 



In order to verify the last statement we first observe 
that 



{0{t)) 



in+tT^t-^ 



?ti 



+ = {^o\RlRte'^^'RlRtORlRte- 



'^+'RlRt\^^o) 
(A24) 



The last equation follows from the definition of time in- 



variance, i^tl^o) = |^o) and RtOR 



t 



O, and from 



the unitarity property R^Rt = 1. Note that equation 



(|A24|) corresponds to the symmetry of time evolutions 
H, which has been discussed previously in 
52l l53l . From the definition of the 7r-boost 



for H ^ - 
References 
we get 



BqU±Bq = -u. 



and can continue equation (A24): 



{0(t))+ = (*o|B, 



h^- 



'^+'BlOB^ 



Qe^«+*B^|*o) 



(A25) 



(0(i))- 



In the last equation we used the assumed 7r-boost invari- 
ance of the observable BqOBq = O and the initial state 
BqI'^o) = l^o), as weh as unitarity Bq = 1. 

In the experiment, we measure the density distribu- 
tion n{rj) = ^cr^jcT^3<y ^^^ ^^^ initial state consists of 
atoms that are completely localized to individual lattice 
sites (see Sec. |A1[ ). Because both the initial state and 
the measured operator are invariant under time-reversal 
and TT-boost, we are guaranteed to find the U ^ —U 
symmetry in the dynamics for all interaction strengths. 

It is interesting to note that the bi-partite character of 
the lattice was crucial to our proof of the symmetry. Thus 
we do not expect to find this symmetry in lattices without 
the bi-partite structure, such as a triangular lattice. 



Non-Interacting atoms: Canceling the harmonic 
confinement 



Figure 12 shows the cloud sizes Rcif) of an expanding 
non-interacting cloud in an 8£^r deep quasi 2D lattice 
as a function of the dipole laser power during the ex- 
pansion. The red line denotes a fit with the expected 
dynamics (Eq. p| to the first 20 ms. While the initial 
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FIG. 12. Expansion of non-interacting atoms as a function of 
dipole beam power p (a.u.) in an %Er lattice. 



expansion velocity depends only slightly on the residual 
confinement, it completely dominates the size after long 
expansion times. 

The largest cloud sizes are reached if the confinement 
created by the dipole trap precisely compensates the an- 
ticonfinement due to the lattice. This situation corre- 
sponds to dipole powers between p = 2/3 and p = 1 in 
Fig. 
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Both an over- and an under-compensation lead 
to deviations from the expected ballistic behaviour and 
limit the cloud size by either classical reflections or Bragg 
reflections of the expanding atoms. 

In the well-compensated case the leading correction to 
the homogeneous Hubbard model (Eq. Al ) arises from 
the fact that the hopping rate in any given direction be- 
comes larger when an atom is not any more in the center 
of the corresponding laser beam. This effect can easily 
be described quantitatively by modeling the experimen- 
tal laser profiles. For a distance of 100 lattice constants 
from the center, the hopping rate increases by about 25%. 
Figure [13] shows the resulting density profiles in the non- 
interacting case, which reproduce all features of Fig. |2] 
of the main text with the exception of the last panel 
at 80 ms, where additional potentials due to technical 
imperfections in the alignment and beam shapes of the 
dipole and lattice beams become important. 

In the interacting case the clouds remain much smaller 
and these effects can completely be neglected. 



5. Image processing and fitting 

In the non-interacting case the perpendicular cloud size 
yj {r'^) is extracted from in-situ phase-contrast images us- 
ing a 2D Gaussian fit: 





t=QhlJ 



t=25hlJ 




t=75hlJ 



FIG. 13. Simulated column density distribution for the non- 
interacting expansion including the finite size of the lattice 
beams. 



Here Xc, ^c, cr^^, a^y, A, 6, are free fit parameters and the 
perpendicular cloud size is given by y/{r'^) = Rg = 
a'^ -\- ay — w'^ ^ where w denotes the imaging resolution 

(radius of Airy disc w < 3/im) of our imaging setup. 

The corresponding single-particle width Rsp{t) is cal- 
culated by deconvolving Rcit) with the initial width 
Rg{0): Rsp{t) = ^Rcitf -RgW' An example of the 
expanded single-particle width is shown in Fig. 3 in the 
main text and agrees well with the linear slope expected 
for a ballistic expansion. 
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FIG. 



14. Numerically calculated density distribution for 
1.2 together with Gaussian fits. 



G(x, y) = Ae 



{y-vcr 



In the case of interacting atoms the shape of the cloud 
changes considerably during the expansion, evolving from 
a "compact" Fermi-Dirac like shape (Supporting Online 
Material of Ref. [6 ) to a "fat tail" distribution, as il- 
lustrated in Fig. [M] using fits to numerically simulated 
data. This leads to considerable systematic errors in the 
estimation of (r^) (7^ R^) in the interacting case. In 
principle, such systematic errors could be avoided by de- 
termining (r^) via direct integration, analogous to the 
(A26) numerical data in Fig. |9] In the experiment, however. 
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this is hindered by imaging aberrations and the smah 
signal to noise ratio in the extreme dilute limit. 

The change in the shape of the cloud is due to the den- 
sity dependent dynamics in the interacting case: While 
the expansion remains ballistic in the low density limit, 
where the mean free path is larger than the distance 
to the cloud edge, the expansion velocity decreases for 
higher densities due to the increasing number of colli- 
sions. As a consequence, (r^) will be dominated by the 
ballistically expanding outermost atoms (see Fig. [6| for 
long expansion times. 
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FIG. 15. Core widths Re as a function of expansion time for 
various interactions in an %Er deep lattice in the 2D case. 
Blue lines den ote th e fits used to extract the core expansion 
velocities (Eq.|A27) 



Rc{t) = ^RM^^vlt^ 



(A27) 



and the resulting core expansion velocities Vc are shown 
in Fig. [5] in the main text. 



6. Width of Feshbach resonance 

The interactions are controlled by use of the well- 
known Feshbach resonance at ^o = 202.1 G between the 
two lowest hyperfine states |F,mF) = |9/2, — 9/2) and 
|9/2, — 7/2) [35 . Compared to previous dipole trap ex- 
periments, the dynamics in the lattice are much more 
sensitive to small scattering lengths, since the strongly 
reduced kinetic energy of atoms in a lattice enhances the 
role of interactions. The observed pronounced depen- 
dence on small interactions (Fig. [5| enabled us to remea- 
sure the zero crossing of the scattering length, which we 
found to be at B{a = 0) = 209.1 ± 0.2 G. Using the 
standard parametrization of the (free space) Feshbach 
resonances 



a{B) = ttbg 1 



w 



B-Bo 

this zero crossing leads to a new width of 
^ = 7.0±0.2G 



(A28) 



(A29) 



compared to the previous dipole trap measurement of 
^dipole = 7.8 ± 0.6 G (PhD. Thesis C. Regal, JILA). 
In addition to the pronounced dependence of the slope, 
only the expansion at the newly assigned zero crossing 
matches that of a single component Fermi gas under the 
same conditions and leads to the square shape expected 
for non-interacting atoms. 



Doublon dissolution time 



In order to focus on the core dynamics we instead use 
the core width Rc^ which denotes the half width at half 
maximum (HWHM) of the azimuthally averaged column 
density distribution, where each distribution is individu- 
ally normalized. This core widths shown in Fig. [15] are 
fitted by the same fit function as in the non-interacting 



The number of atoms on doubly occupied lattice sites is 
measured by converting all doublons into molecules. This 
is achieved by first increasing the lattice depth to 20 Er 
in 200 /is, followed by a magnetic field ramp (5 G/ms) 
over the Feshbach resonance [6 . The resulting doublon 
numbers are fitted by a simple exponential decay. 



